#include <math.h>
#include <R.h>
#include <Rmath.h>
#include <R_ext/Utils.h>
#include "geom3.h"
#include "functable.h"
#include "chunkloop.h"
#include "constants.h"

/*
	$Revision: 1.0 $

	Cross-type pair correlation function of 3D point pattern
	(Epanechnikov kernel) 

	pcf3transCross	  	 translation correction for cross-type
	pcf3isotCross		 isotropic correction for cross-type

  Copyright (C) Adrian Baddeley, Ege Rubak and Rolf Turner 2001-2022
  Licence: GNU Public Licence >= 2
*/

#define FOURPI (2.0 * M_2PI)

void
pcf3transCross(
 Point *p1, int n1,
 Point *p2, int n2,
 Box *b,
 Ftable *pcf,
 double delta
) {
  register int i, j, l, lmin, lmax, maxchunk;
  register double dx, dy, dz, dist;
  register double  vx, vy, vz, tval;
  Point *ip, *jp;
  double dt, vol, lambda1, lambda2, denom;
  double coef, frac, invweight, kernel;

  double sphesfrac(Point *point, Box *box, double r);
  double sphevol(Point *point, Box *box, double r);

  /* compute denominator & initialise numerator*/
  vol = (b->x1 - b->x0) * (b->y1 - b->y0) * (b->z1 - b->z0);
  lambda1 = ((double) n1 )/ vol;
  lambda2 = ((double) n2 )/ vol;
  denom = lambda1 * lambda2;

  for(l = 0; l < pcf->n; l++) {
    (pcf->denom)[l] = denom;
    (pcf->num)[l]   = 0.0;
  }

  /* spacing of argument in result vector */
  dt = (pcf->t1 - pcf->t0)/(pcf->n - 1);

  /* compute numerator */
  OUTERCHUNKLOOP(i, n1, maxchunk, 8196) {
    R_CheckUserInterrupt();
    INNERCHUNKLOOP(i, n1, maxchunk, 8196) {
      ip = p1 + i;
      for(j = 0; j < n2; j++) {
	/* compute pairwise distance */
	jp = p2 + j;
	/* Skip if same point (shouldn't happen for cross-type) */
	if(ip == jp) continue;
	
	dx = jp->x - ip->x;
	dy = jp->y - ip->y;
	dz = jp->z - ip->z;
	dist = sqrt(dx * dx + dy * dy + dz * dz);
	lmin = ceil( ((dist - delta) - pcf->t0) / dt );
	lmax = floor( ((dist + delta) - pcf->t0) / dt );
	if(lmax >= 0 && lmin < pcf->n) {
	  /* kernel centred at 'dist' has nonempty intersection 
	     with specified range of t values */
	  /* compute intersection */
	  if(lmin < 0)
	    lmin = 0;
	  if(lmax >= pcf->n)
	    lmax = pcf->n - 1;
	  /* compute (inverse) edge correction weight */
	  vx = b->x1 - b->x0 - (dx > 0 ? dx : -dx);
	  vy = b->y1 - b->y0 - (dy > 0 ? dy : -dy);
	  vz = b->z1 - b->z0 - (dz > 0 ? dz : -dz);
	  invweight = vx * vy * vz * FOURPI * dist * dist;
	  if(invweight > 0.0) {
	    for(l = lmin; l < pcf->n; l++) {
	      tval = pcf->t0 + l * dt;
	      /* unnormalised Epanechnikov kernel with halfwidth delta */
	      frac = (dist - tval)/delta;
	      kernel = (1 - frac * frac);
	      if(kernel > 0) 	    
		(pcf->num)[l] += kernel / invweight;
	    }
	  }
	}
      }
    }
  }

  /* constant factor in kernel */
  coef = 3.0/(4.0 * delta);

  /* normalise kernel and compute ratio estimate */
  for(l = 0; l < pcf->n; l++) {
    (pcf->num)[l] *= coef;
    (pcf->f)[l] = ((pcf->denom)[l] > 0.0) ?
      (pcf->num)[l] / (pcf->denom)[l] : 0.0;
  }
}

void
pcf3isotCross(
 Point *p1, int n1,
 Point *p2, int n2,
 Box *b,
 Ftable *pcf,
 double delta
) {
  register int i, j, l, lmin, lmax, maxchunk;
  register double dx, dy, dz, dist;
  Point *ip, *jp;
  double dt, vol, denom, mass, tval;
  double coef, frac, kernel;

  double sphesfrac(Point *point, Box *box, double r);
  double sphevol(Point *point, Box *box, double r);

  Point vertex;
  Box   half;

  /* compute denominator & initialise numerator*/
  vol = (b->x1 - b->x0) * (b->y1 - b->y0) * (b->z1 - b->z0);
  denom = ((double) (n1 * n2))/vol;

  for(l = 0; l < pcf->n; l++) {
    (pcf->denom)[l] = denom;
    (pcf->num)[l]   = 0.0;
  }

  /* spacing of argument in result vector */
  dt = (pcf->t1 - pcf->t0)/(pcf->n - 1);

  /* set up for volume correction */
  vertex.x = b->x0;
  vertex.y = b->y0;
  vertex.z = b->z0;
  half.x1  = b->x1;
  half.y1  = b->y1;
  half.z1  = b->z1;
  half.x0  = (b->x0 + b->x1)/2.0;
  half.y0  = (b->y0 + b->y1)/2.0;
  half.z0  = (b->z0 + b->z1)/2.0;

	/* compute numerator */
  OUTERCHUNKLOOP(i, n1, maxchunk, 8196) {
    R_CheckUserInterrupt();
    INNERCHUNKLOOP(i, n1, maxchunk, 8196) {
      ip = p1 + i;
      for(j = 0; j < n2; j++) {
	jp = p2 + j;
	/* Skip if same point (shouldn't happen for cross-type) */
	if(ip == jp) continue;
	
	dx = jp->x - ip->x;
	dy = jp->y - ip->y;
	dz = jp->z - ip->z;
	dist = sqrt(dx * dx + dy * dy + dz * dz);
	lmin = ceil( ((dist - delta) - pcf->t0) / dt );
	lmax = floor( ((dist + delta) - pcf->t0) / dt );
	if(lmax >= 0 && lmin < pcf->n) {
	  /* kernel centred at 'dist' has nonempty intersection 
	     with specified range of t values */
	  /* compute intersection */
	  if(lmin < 0)
	    lmin = 0;
	  if(lmax >= pcf->n)
	    lmax = pcf->n - 1;
	  /* compute edge correction weight for cross-type */
	  /* Use average of the two edge correction weights */
	  mass = 0.5 * (1.0 / sphesfrac(ip, b, dist)
	            + 1.0 / sphesfrac(jp, b, dist));
	  if(mass > 0.0) {
	    mass /= FOURPI * dist * dist;
	    for(l = lmin; l < pcf->n; l++) {
	      tval = pcf->t0 + l * dt;
	      /* unnormalised Epanechnikov kernel with halfwidth delta */
	      frac = (dist - tval)/delta;
	      kernel = (1 - frac * frac);
	      if(kernel > 0) 	    
		(pcf->num)[l] += kernel * mass;
	    }
	  }
	}
      }
    }
  }

  /* constant factor in kernel */
  coef = 3.0/(4.0 * delta);

  /* normalise kernel and compute ratio estimate */
  for(l = 0; l < pcf->n; l++) {
    (pcf->num)[l] *= coef;
    (pcf->f)[l] = ((pcf->denom)[l] > 0.0)?
      (pcf->num)[l] / (pcf->denom)[l]
      : 0.0;
  }
}

/*
  R interface function for cross-type 3D PCF
  This function is called from R via .C()
*/
void
Rcallpcf3cross(
 double *x1, double *y1, double *z1, int *n1,
 double *x2, double *y2, double *z2, int *n2,
 double *x0, double *x1b, double *y0, double *y1b, double *z0, double *z1b,
 double *t0, double *t1, int *m,
 double *f, double *num, double *denom,
 int *method, double *delta
) {
  int i;
  Box b;
  Ftable pcf;
  Point *p1, *p2;
  
  /* Set up box */
  b.x0 = *x0; b.x1 = *x1b;
  b.y0 = *y0; b.y1 = *y1b;
  b.z0 = *z0; b.z1 = *z1b;
  
  /* Set up function table */
  pcf.t0 = *t0;
  pcf.t1 = *t1;
  pcf.n  = *m;
  pcf.f    = f;
  pcf.num  = num;
  pcf.denom = denom;
  
  /* Create point arrays */
  p1 = (Point *) R_alloc(*n1, sizeof(Point));
  p2 = (Point *) R_alloc(*n2, sizeof(Point));
  
  for(i = 0; i < *n1; i++) {
    p1[i].x = x1[i];
    p1[i].y = y1[i];
    p1[i].z = z1[i];
  }
  
  for(i = 0; i < *n2; i++) {
    p2[i].x = x2[i];
    p2[i].y = y2[i];
    p2[i].z = z2[i];
  }
  
  /* Call appropriate function based on method */
  switch(*method) {
    case 0: /* translation correction */
      pcf3transCross(p1, *n1, p2, *n2, &b, &pcf, *delta);
      break;
    case 1: /* isotropic correction */
      pcf3isotCross(p1, *n1, p2, *n2, &b, &pcf, *delta);
      break;
    default:
      error("Unknown correction method");
  }
}